clear all
set more off

global data 	"R:\SharedProjects\Shared2020-070\2016\extend_to_2020\JPE_Replication_dta"
global figures 	"R:\SharedProjects\Shared2020-070\2016\extend_to_2020\JPE_Replication_log"
global in 		"R:\SharedProjects\Shared2020-070\2016\input"



cap log close
log using $figures\P_file,replace t

**************************************************************************************************************
cap program drop pdm1m0
program pdm1m0
	args x1 x2 r
	g pDi=binormal(`x1',`x2',`r')/normal(`x2')						/*This is Pr(D=1|A=1)*/
	g m1i=normalden(`x1')*(normal((1-`r'^2)^(-0.5)*(`x2'-`r'*`x1')))/binormal(`x2',`x1',`r') + ///
		`r'*normalden(`x2')*(normal((1-`r'^2)^(-0.5)*(`x1'-`r'*`x2')))/binormal(`x2',`x1',`r')
	gen m0i=-normalden(`x1')*(normal((1-`r'^2)^(-0.5)*(`x2'-`r'*`x1')))/binormal(`x2',-`x1',-`r') + ///
		`r'*normalden(`x2')*(1-normal((1-`r'^2)^(-0.5)*(`x1'-`r'*`x2')))/binormal(`x2',-`x1',-`r')	
	scalar sigma=sqrt(1+s_csi^2^2)
	gen part1=0.5*(m1i+m0i)/sigma
	gen part2=-sigma*(m1i-m0i)^(-1)*ln(pDi/(1-pDi))
	gen part3=-sigma*(m1i-m0i)^(-1)

	qui su part1 if F==1 & A==1
	scalar p1_1=r(mean)
	qui su part1 if F==0 & A==1
	scalar p1_0=r(mean)
	scalar part1m=p1_1-p1_0
	qui su part2 if F==1 & A==1
	scalar p2_1=r(mean)
	qui su part2 if F==0 & A==1
	scalar p2_0=r(mean)
	scalar part2m=p2_1-p2_0
	qui su part3 if F==1 & A==1
	scalar part3m=r(mean)
	
	qui su pDi if F==1 & A==1	
	scalar pD_1=r(mean)
	qui su pDi if F==0 & A==1	
	scalar pD_0=r(mean)
	qui su m1i if F==1 & A==1	
	scalar m1_1=r(mean)
	qui su m1i if F==0 & A==1	
	scalar m1_0=r(mean)
	qui su m0i if F==1 & A==1	
	scalar m0_1=r(mean)
	qui su m0i if F==0 & A==1	
	scalar m0_0=r(mean)
	
	scalar drop sigma p1_1 p1_0 p2_1 p2_0
	drop pDi m1i m0i part1 part2 part3
end
**************************************************************************************************************

*********************************************************************************************************************************
/*CORRECTION FACTORS*/
cap program drop get_correction
program get_correction
	scalar cfl   =bF_L-((mu-gamma)/(sqrt(1+s_omega^2+s_psi^2))) 
	scalar cfa   =bF_A-((mu-gamma-tau)/(sqrt(1+s_omega^2+s_psi^2+s_v^2)))
	scalar cfr   =bF_R-delta
	scalar cfcla =rho_LA-(sqrt(1+s_omega^2+s_psi^2)/sqrt(1+s_omega^2+s_psi^2+s_v^2))
	scalar cfclr =rho_LR-(-1/(sqrt(1+s_csi^2)*sqrt(1+s_omega^2+s_psi^2)))
	scalar cfcra =rho_RA-(-1/(sqrt(1+s_csi^2)*sqrt(1+s_omega^2+s_psi^2+s_v^2)))
	scalar cfclax=rho_LA/(sqrt(1+s_omega^2+s_psi^2)/sqrt(1+s_omega^2+s_psi^2+s_v^2))
	scalar cfclrx=rho_LR/(-1/(sqrt(1+s_csi^2)*sqrt(1+s_omega^2+s_psi^2)))
	scalar cfcrax=rho_RA/(-1/(sqrt(1+s_csi^2)*sqrt(1+s_omega^2+s_psi^2+s_v^2)))
end
***********************************************************************************************************************************

cap program drop get_baseline_pars
program get_baseline_pars
	preserve
		u $data/R_estimates_DWMD,clear
		destring DWMD,gen(beta)
		replace namev="alpha" if namev=="alpha1"
		ren namev namepar
		*********************************************************************************************************************************
		*BASELINE PARAMETERS
		*********************************************************************************************************************************
		su beta if namepar=="pi"
		scalar mu		=r(mean) 
		su beta if namepar=="gamma"
		scalar gamma	=r(mean)
		su beta if namepar=="tau"
		scalar tau		=r(mean)
		su beta if namepar=="sigma_omega"
		scalar s_omega	=r(mean)
		su beta if namepar=="sigma_psi"
		scalar s_psi	=r(mean)
		su beta if namepar=="sigma_v"
		scalar s_v		=r(mean)
		su beta if namepar=="sigma_csi"
		scalar s_csi	=r(mean)
		su beta if namepar=="alpha"
		scalar alpha1	=r(mean)
		scalar sigma	=sqrt(1+s_csi^2)			/*This is [1+var(noise_SSA)] */	
		***********************************************************************************************************************************
		*********************************************************************************************************************************
	restore
end


cd $data
***********************************************************************************************************************************************************************************
u $data\after_ML, clear

scalar rho_LA=rho_A_L
scalar rho_RA=rho_A_R
scalar rho_LR=rho_L_R
scalar betaL=bF_L
scalar betaA=bF_A
scalar betaR=bF_R

ren female F
gen Aw=1-R
ren logbs Y

ren a1 Lhat										/*This includes the effect of gender*/
ren a2 Ahat										/*This includes the effect of gender*/		
ren a3 Rhat										/*This includes the effect of gender*/

gen xbL=Lhat-bF_L*F			/*subtract the effect of gender because it is going to add it "structurally" later*/
gen xbA=Ahat-bF_A*F			/*subtract the effect of gender because it is going to add it "structurally" later*/
gen xbR=Rhat-bF_R*F			/*subtract the effect of gender because it is going to add it "structurally" later*/	


get_baseline_pars

gen Dhat=xbL*sqrt(1+s_omega^2+s_psi^2)+mu*F
scalar rho_DA=1/sqrt(1+s_omega^2+s_psi^2+s_v^2)

pdm1m0 Dhat Ahat rho_DA

scalar delta=part1m+part2m+part3m*alpha1 	

get_correction

gen cfid=1
replace cfid=sum(cfid)
gen Khat=-Lhat
				
matrix W1=(1, -rho_LA, -rho_LR \ -rho_LA,1,rho_RA \ -rho_LR, rho_RA, 1)
putmata cfid U1=(Khat Ahat Rhat) if !missing(Rhat), replace

mata
	R1=vech(st_matrix("W1"))'
	Z=mvnormal(U1,R1)
	st_matrix("Z",Z)
end

getmata Z, id(cfid) replace
g Type2_baseline=1-(Z/binormal(Ahat,Khat,-rho_LA))
su Type2* if F==1 & A==1
su Type2* if F==0 & A==1
save trial_type2,replace


******pi*******
get_baseline_pars

matrix define mu_results=(0,0,0,0,0,0,0,0,0)

forvalues j=-1(0.1)1	{
	u trial_figure1,clear
	scalar mu_x=`j'

	g Lhat_pi=xbL+(((mu_x-gamma)/(sqrt(1+s_omega^2+s_psi^2)))+cfl)*F
	g Ahat_pi=xbA+(((mu_x-gamma-tau)/(sqrt(1+s_omega^2+s_psi^2+s_v^2)))+cfa)*F
	g Dhat_pi=/*xbLV*/xbL*sqrt(1+s_omega^2+s_psi^2)+mu_x*F
	scalar rho_DA=1/sqrt(1+s_omega^2+s_psi^2+s_v^2)
	pdm1m0 Dhat_pi Ahat_pi rho_DA
	scalar delta=part1m+part2m+part3m*alpha1
	g Rhat_pi=xbR+(delta+cfr)*F
	gen Khat_pi=-Lhat_pi

	matrix W1=(1, -rho_LA, -rho_LR \ -rho_LA,1,rho_RA \ -rho_LR, rho_RA, 1)			
	putmata cfid U1=(Khat_pi Ahat_pi Rhat_pi) if !missing(Rhat_pi), replace

	local endmata end
	mata
		R1=vech(st_matrix("W1"))'
		Z=mvnormal(U1,R1)
		st_matrix("Z",Z)
	`endmata'
	getmata Z, id(cfid) replace
	
	g Type2_pi=1-(Z/binormal(Ahat_pi,Khat_pi,-rho_LA))

	qui su Type2_pi if F==1 & A==1			
	scalar t1f=r(mean)
	qui su Type2_pi if F==0 & A==1
	scalar t1m=r(mean)
	*********************************************************************************************************************************
	matrix mu_results=mu_results\(`j',t1f,t1m,m1_1,m1_0,m0_1,m0_0,pD_1,pD_0)
	*********************************************************************************************************************************
	}

svmat mu_results
keep if mu_results5!=.
keep mu_result*
drop if mu_results5==0
ren mu_results2 	pi_Type2_F
ren mu_results3		pi_Type2_M
ren mu_results4		pi_m1_F
ren mu_results5		pi_m1_M
ren mu_results6		pi_m0_F
ren mu_results7		pi_m0_M
ren mu_results8		pi_pD_F
ren mu_results9		pi_pD_M
ren mu_results1 	value
sort value
save mu_results_type2,replace


******gamma*******
get_baseline_pars

matrix define gamma_results=(0,0,0,0,0,0,0,0,0)

forvalues j=-1(0.1)1	{
	u trial_figure1,clear
	scalar gamma_x=`j'

	g Lhat_gamma=xbL+(((mu-gamma_x)/(sqrt(1+s_omega^2+s_psi^2)))+cfl)*F
	g Ahat_gamma=xbA+(((mu-gamma_x-tau)/(sqrt(1+s_omega^2+s_psi^2+s_v^2)))+cfa)*F
	g Dhat_gamma=xbL*sqrt(1+s_omega^2+s_psi^2)+mu*F
	scalar rho_DA=1/sqrt(1+s_omega^2+s_psi^2+s_v^2)
	pdm1m0 Dhat_gamma Ahat_gamma rho_DA
	scalar delta=part1m+part2m+part3m*alpha1
	g Rhat_gamma=xbR+(delta+cfr)*F
	gen Khat_gamma=-Lhat_gamma

	matrix W1=(1, -rho_LA, -rho_LR \ -rho_LA,1,rho_RA \ -rho_LR, rho_RA, 1)			
	putmata cfid U1=(Khat_gamma Ahat_gamma Rhat_gamma) if !missing(Rhat_gamma), replace

	local endmata end
	mata
		R1=vech(st_matrix("W1"))'
		Z=mvnormal(U1,R1)
		st_matrix("Z",Z)
	`endmata'
	getmata Z, id(cfid) replace
	
	g Type2_gamma=1-(Z/binormal(Ahat_gamma,Khat_gamma,-rho_LA))

	qui su Type2_gamma if F==1 & A==1			
	scalar t1f=r(mean)
	qui su Type2_gamma if F==0 & A==1
	scalar t1m=r(mean)
	*********************************************************************************************************************************
	matrix gamma_results=gamma_results\(`j',t1f,t1m,m1_1,m1_0,m0_1,m0_0,pD_1,pD_0)
	*********************************************************************************************************************************
	}

svmat 	gamma_results
keep if gamma_results5!=.
keep 	gamma_result*
drop if gamma_results5==0
ren gamma_results2 		gamma_Type2_F
ren gamma_results3		gamma_Type2_M
ren gamma_results4		gamma_m1_F
ren gamma_results5		gamma_m1_M
ren gamma_results6		gamma_m0_F
ren gamma_results7		gamma_m0_M
ren gamma_results8		gamma_pD_F
ren gamma_results9		gamma_pD_M
ren gamma_results1 	value
sort value
save gamma_results_type2,replace


******tau*******
get_baseline_pars

matrix define tau_results=(0,0,0,0,0,0,0,0,0)

forvalues j=-1(0.1)1	{
	u trial_figure1,clear
	scalar tau_x=`j'
	g Lhat_tau=xbL+(((mu-gamma)/(sqrt(1+s_omega^2+s_psi^2)))+cfl)*F
	g Ahat_tau=xbA+(((mu-gamma-tau_x)/(sqrt(1+s_omega^2+s_psi^2+s_v^2)))+cfa)*F
	g Dhat_tau=xbL*sqrt(1+s_omega^2+s_psi^2)+mu*F
	scalar rho_DA=1/sqrt(1+s_omega^2+s_psi^2+s_v^2)
	pdm1m0 Dhat_tau Ahat_tau rho_DA
	scalar delta=part1m+part2m+part3m*alpha1
	g Rhat_tau=xbR+(delta+cfr)*F
	gen Khat_tau=-Lhat_tau

	matrix W1=(1, -rho_LA, -rho_LR \ -rho_LA,1,rho_RA \ -rho_LR, rho_RA, 1)			
	putmata cfid U1=(Khat_tau Ahat_tau Rhat_tau) if !missing(Rhat_tau), replace

	local endmata end
	mata
		R1=vech(st_matrix("W1"))'
		Z=mvnormal(U1,R1)
		st_matrix("Z",Z)
	`endmata'
	getmata Z, id(cfid) replace
	
	g Type2_tau=1-(Z/binormal(Ahat_tau,Khat_tau,-rho_LA))

	qui su Type2_tau if F==1 & A==1			
	scalar t1f=r(mean)
	qui su Type2_tau if F==0 & A==1
	scalar t1m=r(mean)
	*********************************************************************************************************************************
	matrix tau_results=tau_results\(`j',t1f,t1m,m1_1,m1_0,m0_1,m0_0,pD_1,pD_0)
	*********************************************************************************************************************************
	}

svmat 	tau_results
keep if tau_results5!=.
keep 	tau_result*
drop if tau_results5==0
ren tau_results2 		tau_Type2_F
ren tau_results3		tau_Type2_M
ren tau_results4		tau_m1_F
ren tau_results5		tau_m1_M
ren tau_results6		tau_m0_F
ren tau_results7		tau_m0_M
ren tau_results8		tau_pD_F
ren tau_results9		tau_pD_M
ren tau_results1 	value
sort value
save tau_results_type2,replace


******alpha*******
get_baseline_pars

matrix define alpha_results=(0,0,0,0,0,0,0,0,0)

forvalues j=-1(0.1)1	{
	u trial_figure1,clear
	scalar alpha_x=`j'
	g Lhat_alpha=xbL+(((mu-gamma)/(sqrt(1+s_omega^2+s_psi^2)))+cfl)*F
	g Ahat_alpha=xbA+(((mu-gamma-tau)/(sqrt(1+s_omega^2+s_psi^2+s_v^2)))+cfa)*F
	g Dhat_alpha=xbL*sqrt(1+s_omega^2+s_psi^2)+mu*F
	scalar rho_DA=1/sqrt(1+s_omega^2+s_psi^2+s_v^2)
	pdm1m0 Dhat_alpha Ahat_alpha rho_DA
	scalar delta=part1m+part2m+part3m*alpha_x
	g Rhat_alpha=xbR+(delta+cfr)*F
	gen Khat_alpha=-Lhat_alpha
	matrix W1=(1, -rho_LA, -rho_LR \ -rho_LA,1,rho_RA \ -rho_LR, rho_RA, 1)			
	putmata cfid U1=(Khat_alpha Ahat_alpha Rhat_alpha) if !missing(Rhat_alpha), replace

	local endmata end
	mata
		R1=vech(st_matrix("W1"))'
		Z=mvnormal(U1,R1)
		st_matrix("Z",Z)
	`endmata'
	getmata Z, id(cfid) replace
	
	g Type2_alpha=1-(Z/binormal(Ahat_alpha,Khat_alpha,-rho_LA))

	qui su Type2_alpha if F==1 & A==1			
	scalar t1f=r(mean)
	qui su Type2_alpha if F==0 & A==1
	scalar t1m=r(mean)
	*********************************************************************************************************************************
	matrix alpha_results=alpha_results\(`j',t1f,t1m,m1_1,m1_0,m0_1,m0_0,pD_1,pD_0)
	*********************************************************************************************************************************
	}

svmat 	alpha_results
keep if alpha_results5!=.
keep 	alpha_result*
drop if alpha_results5==0
ren alpha_results2 		alpha_Type2_F
ren alpha_results3		alpha_Type2_M
ren alpha_results4		alpha_m1_F
ren alpha_results5		alpha_m1_M
ren alpha_results6		alpha_m0_F
ren alpha_results7		alpha_m0_M
ren alpha_results8		alpha_pD_F
ren alpha_results9		alpha_pD_M
ren alpha_results1 	value
sort value
save alpha_results_type2,replace




u mu_results_type2,clear
merge value using gamma_results_type2
drop _merge
sort value
merge value using tau_results_type2
drop _merge
sort value
merge value using alpha_results_type2
drop _merge
sort value


foreach x in Type2_F Type2_M m1_F m1_M m0_F m0_M pD_F pD_M {
	lab var pi_`x' "{&pi}"
	lab var gamma_`x' "{&gamma}"
	lab var tau_`x' "{&tau}"
	lab var alpha_`x' "{&alpha}"
	}
lab var value "Parameter value"

scatter *_Type2_F value,c(l l l l) s(i o dh +) clp(dash solid solid solid) legend(col(4) pos(6)) graphr(c(white)) ytitle(Pr(Type II error)) ///
	xlabel(-1(.2)1) 
graph export $figures\figureOA3.pdf,replace
	
erase trial_figure1.dta
erase trial_type2.dta

